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As one approaches the continuum limit, QCD systems, investigated via numerical 
simulations, remain trapped in sectors of field space with fixed topological charge. As 
a consequence the numerical studies of physical quantities may give biased results. 
The same is true in the case of two dimensional CP^~^ models. In this paper we 
show that metadynamies, when used to simulate CP^~^, allows to address efficiently 
this problem. By studying CP^^ we show that we are able to reconstruct the free 
energy of the topological charge F{Q) and compute the topological susceptibility as 
a function of the coupling and of the volume. This is a very important physical 
quantity in studies of the dynamics of the 9 vacuum and of the axion. This method 
can in principle be extended to QCD applications. 
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I. INTRODUCTION 


In numerical simulations of asymptotically free theories, CP^~^ in two dimensions and 
SU{N) in four dimensions, one observes an increase of the autocorrelation time r, defined 
as the number of iterations needed to generate independent field configurations, as one 
proceeds towards the continuum limit where the coupling constant vanishes and the typical 
correlation length diverges. This corresponds to the critical slowing down occurring in 
statistical systems close to a second order phase transition. In addition, in these systems, 
a particularly dramatic increase of autocorrelation time is observed in the case of the 
topological charge, independently of the precise discretized definition which is used in the 
simulation on the lattice [iHin]. 

In general, the asymptotic scaling behavior of the autocorrelation time with the lattice 
spacing is expected to be power-like for all the quantities. On the contrary, in both CP^~^ 
and QCD the autocorrelation time of the topological charge, tq, grows so rapidly with the 
length scale that an apparent exponential behavior tq ~ exp(c^®) in the explored range 
of values of ^ is compatible with the available data W- On the other hand, this peculiar 
effect is not observed for “quasi-Gaussian modes” such as the plaquette or the Polyakov 
line correlations, suggesting a separation of the dynamics of the topological modes from 
that of quasi-Gaussian ones. 

The different dynamical behavior of quasi-Gaussian and topological modes is induced by 
sizable free-energy barriers separating different regions of the configuration space. There¬ 
fore the evolution in this space presents a long relaxation time due to the transitions 
between different topological charge sectors, and the corresponding autocorrelation time is 
expected to behave as tq ~ exp(AF(Q)), where AF{Q) is the typical free-energy barrier 
between different topological sectors. If the height of the barrier increases as we proceed 
toward the continuum limit, the system can be trapped in a topological-charge sector for 
a number of simulation steps comparable or even longer than the total available resources 
of the simulation. In QCD this problem may bias the lattice predictions for several impor¬ 
tant physical quantities such as the mass of the rj' meson [TTIIIT] . or the polarized baryon 
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structure function [T5I - I18) . In general this bias will be present for any observable correlated 
to the topological charge, whenever the algorithm does not explore with the appropriate 
weight the different topological sectors. This has been shown to occur in ref. [9] in the 
contest of CP^~^ models. It is clear that in order to design a cure for these problems one 
should also measure and understand how the typical free-energy barriers scale with the 
correlation length (and with the physical volume) in a simulation. 

In the recent literature the problem of very long autocorrelation times and of the theo¬ 
retical control over the systematic error in numerical simulations has been addressed in a 
series of papers PElElIMn] and several solutions to the topological critical slowing down 
have been proposed H El II9H25]. 

In this paper we propose to address the problem of topological trapping by using meta¬ 
dynamics, a powerful method that was introduced to enhance the probability of observing 
rare conformational changes and reconstructing the free energy in biophysics, chemistry 
and material sciences systems I261E7]. 

The metadynamics method I261E7] requires the preliminary identihcation of Collective 
Variables (CV) which are assumed to be able to describe the process of interest. A CV 
is a (physical) quantity which depends on a set of the coordinates or fields of the system. 
In the simplest case the variable could be just the topological charge itself although the 
full power of this approach relies in its ability to treat several CVs simultaneously. The 
dynamics in the space of the chosen CVs is enhanced by a history-dependent potential 
constructed as a sum of Gaussians centered along the trajectory followed by the CVs. The 
sum of Gaussians is exploited to reconstruct iteratively an estimator of the free energy. 
This procedure resembles the Wang and Landau algorithm |28], in which a time-dependent 
bias is introduced to modify the density of states to produce flat histograms in models with 
discrete energy levels. 

The history of microscopic systems in normal Monte Garlo or Hybrid Monte Garlo 
(HMG) [To] corresponds to a random walk with a bias in the direction of lower free energy. 
In systems with many local minima the probability to explore the transition regions and 
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tunnel in a different minimum is very small, particularly when the height of the barrier 
is high. In metadynamics, the system has access to a feedback which during the time 
evolution hlls the local free energy minima. Thus, even if at the beginning the system 
visits more often the region at the bottom of a local minimum, after a few steps almost 
deterministically it starts exploring regions corresponding to higher and higher values of 
the free energy. Sooner or later, the system hlls the minimum, climbs out of it, and visits 
another minimum that is eventually also hlled, until all the relevant minima are explored. 
The key idea of metadynamics is exploiting the time-dependent bias potential itself as a 
free energy estimator. In particular, the time average of the bias potential has been shown 
to converge to the negative of F with an error that scales to zero with the inverse square 
root of the simulation time |29) . 


We here propose to address the problem of topological trapping by performing meta¬ 
dynamics on the topological charge. We will show that this approach induces a large 
number of transitions between different sectors, and therefore converges very rapidly. At 
the same time, the approach allows computing the unbiased average value of any observ¬ 
able by standard reweighting techniques. In order to test our proposal we hrst study the 
two-dimensional CP^~^ models that have several features in common with QCD, such 
as asymptotic freedom and a non-trivial topological structure. Since these models require 
much smaller computing resources, they are an ideal theoretical laboratory to be used in 
an early and exploratory stage of any new algorithm. We hnd the improvement induced by 
metadynamics considerable and worth to be implemented in a QCD study that we plan 
to perform in the near future. 


The paper is organized as follows: in Sec. [^we recall the basic ingredients of CP 


N-l 


and define the quantities that will be measured in our numerical simulation; in Sec. Ill 


the metadynamics algorithm is described; in Sec. IV we present the results of our numer¬ 


ical study and compare several quantities obtained without or with metadynamics. Our 
conclusions will be given in Sec. |V] and some technical details are presented in appendix [A| 




II. CP’^-^ AND THE DIFFERENT DEFINITIONS OF THE TOPOLOGICAL 

CHARGE ON THE LATTICE 
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In the continuum the two dimensional CP^ ^ model is dehned by the action 

d'^xD^zDf.z, ( 1 ) 

where z s, complex iV-dimensional held with unit norm ^ ~ 1 

covariant derivative is given by = d^ + iA^. The held has no kinetic term. For this 
reason, by using the equation of motion, can be expressed in terms of the held 2 : 

^ (f ■ d^z - d^z ■ z) . (2) 

The action is invariant under the local U{1) gauge transformation 

z{x) pAP ^ Af,- d^A{x) , ( 3 ) 

The connected correlation function and magnetic susceptibility are dehned as 

G{x) = {TrP{x)P{0)) — ^(rrP(a;))(rrP(0)), P{x) = z{x) z{x) {Pij{x) = z*Zj) 

Am = J d‘^xG{x) = G{p = 0). (4) 

In the continuum one can dehne the correlation length ^ from the large time-distance 
behavior of G{x) 

G{xo) = j dxiG{xQ,Xi) ~ e~^o/iw ^ (' 5 ^ 

or from the second moment of the correlation function 

fd^xfG(x) 
f d?xG{x) 


( 6 ) 
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In our numerical analysis we have used the dehnition given below in eq. (12), which is 
proportional to and in the scaling region. 


The topological-charge density g(x), the total topological charge Q and the topological 
susceptibility Xt are dehned as 

q{x) = Q = J d‘^xq{x), 

Xt = J (fx{q{x)q{0)) = ■ (7) 

The lattice action is given by [30] 

S = - ^ D^zhD^zh , ( 8 ) 

d n,fj, 

where we introduce the lattice covariant derivative 




(9) 


expressed in term of the U{1) gauge link = exp{iA^{n -|- fta/2), where a is the lattice 
spacing and jl the unit vector in the direction /i. This is the simplest nearest neighbor 
derivative. Improved versions of it could also be used. Lattice gauge invariance, in its 
linear realization, reads 


~ ^ 1 '^n ~ ^ 


( 10 ) 


By integrating by part the term in eq. ([^ and removing the term independent of the helds 
we arrive to the action used in our numerical simulation 


* 5 * ^ Zn+fi^n,ji^n 1 ] • 


( 11 ) 


We use the lattice dehnition of the correlation length ^g [2] 


— A ^-2 


G(0) - G(g, 


4 sin {q^/2) G'(g^) 


(12) 
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where qm is the smallest non zero dimensionless momentum on a lattice with lattice spacing 
a and physical volume aL, namely qm = (27r/L,0). In the continuum limit a^c iw) 
where the correlation was dehned in eq. ([^. 

At hnite lattice spacing no unique dehnition of the topological charge exists. One 
possible geometrical dehnition of the topological charge was introduced in Ref. j3T] as 


= — y^ImjlnTr [P(n + /t + i>)P(n + /i)P(n)] + liiTr [P(n + i>)P(n + fi + z>)P(n)]} , 

271 

n 

(13) 

where ^ ^ v and P is the projector dehned in eq. Q. Thanks to the periodicity of 
the lattice, the previous dehnition is guaranteed to be an exact integer number on every 
conhguration (strictly speaking this in fact holds for all conhgurations except for a subset 
of measure zero). Diherent geometrical dehnitions of the topological charge, however, are 
not guaranteed to return the same integer number, especially on coarse lattices. can 
be written in terms of the phase 6*^^^ = arg {zfiZfi+f,) as 

^ ^ {On,fL + dft+fifi - On+i>,il - On,u) ( 14 ) 


Being strictly integer-valued, the geometrical dehnition of eq. (13) is insensible to in- 
hnitesimal continuum deformations of the helds. As a side ehect, a geometrical version 
of the topological charge cannot be used to dehne a bias variable in the metadynamics 
approach. Indeed the HMC algorithm dehnes a continuous dynamics in a hctitious-time, 
which requires the evaluation of the derivative of the action with respect to the local vari¬ 


ables. Any bias potential built in terms of the geometrical dehnition of eq. (13) would 
give rise to a dynamics completely insensible to the past history of the topological charge, 
as long as the system remains conhned in the same topological sector. The system would 
only feel the ehects of the bias potential when crossing the boundaries between diherent 
topological sectors, where the inhnite force arising from the discontinuity in the potential 
would break the numerical integration of the equations of motion. 




8 


Another definition of the topological charge that better serves onr scope is given in 
terms of the imaginary part of the plaquette, as illustrated in Ref. (3^ 

+ fi)X^{n +i))X^{n)] < u. (15) 

n 

differs from the geometrical definition of eq. ( [T^ . We have 

Q^ = ZQQ^ + rj, (16) 

where Zq is a suitable renormalization constant that can be computed in perturbation 
theory and r] is a zero-average additive ultraviolet noise depending on the particular field 
configuration, the variance of which is a increasing function of the lattice volume and a 
mildly dependent function of the lattice spacing. 


Such noise can be reduced by computing the topological charge of eq. (15) after “smooth¬ 


ing” the A fields (see e.g. Ref. for a recent comparison between cooling and Wilson 
flow in the contest of QCD), or by smearing the gauge fields with procedures like APE j35] 
or HYP m. 


Introducing a bias potential related to a non-geometric definition of the topological 
charge may accelerate the dynamics of the ultraviolet noise, that is expected to be con¬ 
nected with the degrees of freedom ultimately related to the tunneling between different 
topological sectors. As a matter of fact, it will turn out to be useful to filter out the high¬ 
est frequency components of the noise, with the purpose of capturing those components 
expected to be more closely related to the tunnelling phenomenon. To this end we make 
use of a modification of the stout smearing 1371 , adapted to treat U{1) variables, according 
to the procedure explicitly described in appendix 

We stress that the smeared topological charge is exploited in our approach only to 
enhance the tunneling between different minima, thus improving the convergence of the 
system to thermalisation with the possibility of averaging over all the different topological 
sectors. After convergence is achieved, one can then compute the properties of the system 
described by the normal action by using standard reweighting techniques. For 

example, one can compute the probability of observing different values of the exact Q. 
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To conclude this section we define the fictitious-time autocorrelation function Co{t), 
where O is any observable (topological charge, magnetization, ^ etc.) and t is the discrete 
simulation time expressed in sweeps 

Co(t) = ((O(t)-(O))(O(0)-(O))), (17) 

where the averages are taken after the thermalisation of the system. We expect that 
Coit) e at large f, where To is the typical autocorrelation time of the observable 

O. 


III. METADYNAMICS 


Let us consider a physical system described by a set of coordinates x and an action S{x) 
that evolves under the effect of a given dynamics that samples a equilibrium probability 
proportional to exp(—5'(a;)). We are interested in exploring the properties of the system 
as a function of a CV Q{x). The probability distribution of this variable is given by 

P{q) = J dx exp{—S{x)) 5{q — Q{x)). (18) 

and the corresponding free energy is F [q) = — log(P(g)). The probability distribution 
and the free energy can be estimated by computing the histogram of g = Q{x) over a phase 
space trajectory x{t): P{q) rs-/ \ Jq dt' 5{q {t') — q) where q{t) = Q {x (t)). If the system 
displays metastability, namely if P{q) is characterized by the presence of two or more local 
maxima, separated by regions where the probability of observing the system is small, this 
estimator is affected by large statistical errors, since q is normally trapped in a region of 
high probability, and only rarely performs a transition to another region. This is the case 
for the dynamics of a CP^~^ model, where at small lattice spacing the system is trapped 
in a specific topological sector. 

In metadynamics the action S (x) is modified by adding to it a history-dependent po- 
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tential of the form 

VG{Q{x),t) = 9 {Q (x) - q {t')) (19) 

t' = Tg,2tg,.. . 
t' < t 

where g (g) is a non-negative function of its argument, that rapidly vanishes for large |g|. 
In the original implementation, g{q) = wexp ''^here w and 6q are two parameters 

that can be used to tune the growth speed of Vg- Thus, the metadynamics potential is a sum 
of small repulsive potentials, acting on the CV, and localized on all the conhgurations q{t) 
that have already been explored during the trajectory, up to time t. This potential disfavors 
the system from revisiting conhgurations that have already been visited. If the dynamics 
of q is bound in a hnite connected region, after a transient the probability distribution 
of q in this region can only be approximately hat. Indeed, if this is not the case, by 
dehnition the system would spend more time in a subregion of q, and Vg would grow 
preferentially in that region, disfavoring the system from remaining there. Thus, deviations 
from the hat distribution can only survive for a transient time. P{q) exp {—Vg {q, t)) must 
be approximately constant or, equivalently, 

ydqV) - P{q)- (20) 


This equation states that in metadynamics the free energy is estimated by the negative 


of the bias potential itself. More precisely, since eq. (20) is valid at any time, the best 
estimator of the free energy at time t is given by the (large) time average of Vg up to time 

t, 

_ 1 P 

( 21 ) 


F(q) ~ Va(q,t) = 


t — t, 


dt'VG{q,t') 


eq Jte 


The equilibration time teq entering in eq. (21) is the time at which the history dependent 
potential becomes approximately stationary (or, equivalently, the probability distribution 
as a function of q becomes approximately hat). Like in the ordinary estimates of the average 
value of an observable, the exact choice of t^g inhuences only the convergence speed, but 


not the hnal result. The diherence between —F and Vg in eq. (21) decreases as the square 
root of f — teq, with a prefactor that strongly depends on the specihc CV q m- 
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Since the bias potential in eq. (19) must be a differentiable function of the fields, we 


use as CV the topological charge defined in eq. (15) rather than the “integer" charge of 


eq. (13) 


The probability of observing different values of the exact Q is then computed by 
reweighting, as discussed below. 


The most efficient choice of the smoothing parameters is such that the distribution 
probability of is the largest which still allows an assignment to different integer values 
of the topological charge. In other words, the distribution probabilities of for different 
given topological sectors of are such that they do not overlap. This determines the 
optimal width of Q^. 


For the sake of computational efficiency, we store the history-dependent potential on a 
regular grid of spacing 6q; (go, <?i, • • ■ , <?n), with g* = go -|- i 6q. The use of the grid makes it 
possible to carry on metadynamics for long runs at a fixed overhead per sweep (in computer 
time), whereas the computer time with the naive procedure would linearly increase with 
the number of sweeps. 

At the beginning of the simulation we set Vi = Vg (g*) = 0. Then, at every step, we 


1. compute the value of the CV q{t) = Q^(t); 

2. find the grid interval i where it falls 


i = int 


f q ji) - % 

V h 


+ 0.5 


3. update the history dependent potential as follows 

q (t) - qi 


V = V + n; 1- 


Sq 


Vi+i — V)+i + w 


q ji) - qi 

6q 


( 22 ) 


Thus, in our implementation the function g{q) entering in eq. (18) is a sum of two 
triangular-shaped functions, one centered on g*, the other on g^+i. 
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The force ruling the evolution of the helds x is then changed by adding to it the com¬ 
ponent deriving from the history-dependent potential. Thus, the total force is 


VL (x) - VVG{Q{x),t) = -VL {x) - + (x) 

6q 


(23) 


The two crucial parameters in this procedure are w and 6q. For too large values of 
w an accurate integration of the trajectory would require an inhnitesimal time step (for 
vanishing w metadynamics becomes the standard HMC). The optimal grid spacing Sq must 
be such that i) the potential wells are hlled rapidly, and this requires a large 6q] ii) the free 


energy F{q), eqs. (18) and (20), can be accurately reconstructed, and this requires a small 


value of 6q. Suitable values for w and 6q can be initially estimated in an unbiased run by 
measuring the transition probability between two topological sectors and the fluctuations 
of in a given sector respectively [27]. The hrst will give us the height of the barrier, 
and w must be significantly smaller than this height, while second give us an estimate of 
the width of the distribution. 

In order to reach a stationary state in which the probability distribution as a function of 
q is flat it is necessary constraining the dynamics of g in a hnite region. This can be done 
without loss of generality by suppressing conhgurations, corresponding to large values of q, 
that have a exponentially low probability weight. A well established method iza is given 
by the introduction of a threshold value of the topological charge we disfavor values 
of larger than the threshold by a repulsive force which increases linearly starting from 

VrestiQix)) = k {Q{x) - , \Q{x)\ > Q*"'’. (24) 

A; is a new parameter which should be chosen of the same size as w/{25q^). can be 
estimated from the value of the topological susceptibility, choosing Qthr 3> ^XqL'^- 

In summary hve parameters have to be tuned, namely w, Sq, , k and the smearing 


parameter pn discussed in appendix A. The grid on q in eq. 22 is chosen in such a way that 
qi = —and g„_i = The forces from metadynamics in eq. 23 are added only if 


Q{x) E [—-|-(5*"’'].The history-dependent potential Vi is updated according to eq. 22 
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also when Q(x) G [go, gi] or Q(x) G [gn-i, Qn]- Note that, since the threshold is chosen 
in such a way that the conhgurations that are disfavored by the restraining potential in 
eg. [^correspond to a exponentially low probability, using a larger value of would not 
change the estimate of the average value of any observable, except for exponentially small 
corrections, much smaller than the statistical error. 

More details on how the parameters of this approach are chosen and tuned to optimize 
the efficiency can be found in ref m 


IV. NUMERICAL RESULTS 


In this section we present a comparison of results obtained by using the standard HMC 
algorithm with those obtained by using metadynamics. In particular we focus on the 
autocorrelation time r of several physical quantities among which the autocorrelation of 
the topological charge tq and its scaling properties in the continuum limit. 

We have studied with N = 21, at several values of the coupling constant, 

with different physical volumes at fixed correlation length and with different values of the 
correlation length at fixed physical size. We have chosen because it was already 

extensively studied in the literature, see for example refs. PEI [32]. These papers made 
the choice of a large value of N in order to enhance hnite size effects and the critical slowing 
down that they wanted to study. We have taken the same value of N for the same reasons 
and in order to have a direct comparison with previous, quite precise data regarding the 
main observables (energy, correlation length, magnetic and topological susceptibility). 

We performed molecular dynamics with a time-step At = 1 in the hctitious-time, in¬ 
tegrating the equation of motion by means of the Omelyan integrator, using 18 steps per 
trajectory to achieve a high acceptance rate (0(90%)). 

We present in Tables |A| and |B| the input parameters (coupling /?, number of sweeps W, 
lattice size L) and the results of the numerical simulations for the following observables: 
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the dimensionless correlation length dehned in eq. ( [1^ and L/^g] the average energy 
density E = gS/2V with the action S dehned in eq. 0; the magnetic susceptibility Xm 
and the topological susceptibility xq- We have chosen the lattice volume in order to work 
with hxed physical hnite volumes at the different values of f3 (^g); namely at fixed L/^g- 
For comparison, in tables and we present the values for ^Gi Xm and Xt obtained with 
the standard HMC and with metadynamics, respectively. As for the relevant parameters 


of the metadynamics, see eq. (22), we used 5q = 0.05 and varied tn as a function of the 


coupling {w = 0.025 at /3 = 0.65 and w = 0.140 ai (3 = 0.65) to compensate the increase 
of the height of the barriers in the continuum limit. 

The unbiased expectation value of an observable O is computed through the reweighting 
procedure 


( 0 ) = 


E.O.exp(-y(Q>)) 

E.exp(-V'(Q<')) 


(25) 


where V{Q^)) is defined in eq. (21). 


Since this work describes the first application of metadynamics to lattice field theory, 
in the next subsection we will present several details about the numerical results, whereas 


in Sec. (IVB) we discuss the scaling properties of the efficiency as we proceed toward the 


continuum limit. 


A. A comparison of standard HMC and metadynamics 

We start by considering a metadynamics run performed at /3 = 0.70 and = L = 62. For 
these values, the standard HMC is still able to explore different topological sectors, and 
achieve convergence. This allows checking if the two approaches provide consistent results. 

In Fig. [^we show the metadynamics bias potential averaged over progressively larger 
number of MC sweeps. As the simulation proceeds, the minima are iteratively filled, until 
all the five minima shown in figure are explored, and the bias potential starts growing evenly 
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FIG. 1. Running average of the metadynamics bias potential VciQjt) defined in eq. (19) for a 
run performed with fi = 0.70 and L — 62. The average is taken over 1000 sweeps, and shown at 
intervals of 1000 sweeps, starting from the bottom to the top. 


in all the Q range. At this point, the free energy estimator Vc{q,t) can be considered 
converged. Since the instantaneous value of Vci^q^t) is subject to fluctuations, in order 
to extract the free energy and its statistical error, we average Vciq,!) after convergence, 
namely after 10k HMC sweeps for the example in Figure. In Fig. the reconstructed 
average free energy of the topological charge F (Q^) and its statistical uncertainty, shown 
as a (orange) band, estimated by dividing the Monte Carlo history after equilibration 
in Hp = 4 different intervals and estimating the standard deviation of the np measures. 
In the bottom panel we compare the metadynamics result with —log{P{Q^)) estimated 
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in a standard HMC run. Remarkably, the two estimates are fully consistent within the 
small error bars, indicating that metadynamics allows computing reliably the probability 
distribution of the charge. 

The most important effect of the metadynamics bias is reducing the autocorrelation 
time of the observables by orders of magnitude. In Fig. we show the (hctitious) time 
autocorrelation function for the magnetization and topological charge for the standard 
HMC run and for the metadynamics run. For comparison, we have chosen for these hgures 
a value of f3 and of the lattice size L equal to one of the runs presented in Ref. |32] 
(see also m)- The improvement is signihcant: by htting the exponential decrease of the 
autocorrelation functions we hnd tq = 155000 ± 10000 in the HMC run, whereas the 
corresponding value for the metadynamics run is tq = 5600 ± 1000. In both cases we 
hnd Tm = 50 ± 10. A related quantity is the transition probability per unit time between 
two different topological sectors, u. This quantity is dehned as the number of jumps 
between two different values of divided by the total number of sweeps. For the two 
runs corresponding to Fig. 1^ we have u = 4.98 ■ 10 ® with the HMC and u = 2.24 ■ 10 ^ 
for metadynamics respectively. 

One of the main results of this paper is the reconstructed free energy. Fig. at different 
values of the coupling (3 and of the volume. Although the specihc form of the free energy 
depends on the dehnition of Q^, it contains some important physical information. F{Q) is 
very well approximated by the function 

F(Q) = + 5 sin^(7rQ), (26) 

where A and b are numerical constants which in general, at hxed N, depend on L and (3. 

As an example we htted the effective potential reported in the top-left panel of Fig. 
relative to f3 = 0.70 with L = 62. Taking c = 0, the results of the £t, in the metadynamics 
case, are A = 0.20 ± 0.06 and B = 4.38 ± 0.10, with a x^/dof 1. By using the relation 
A ~ l/(2xtR), this corresponds to Xt = (6.50 ± 1.9) x 10“^ well compatible with the 
results given in Tab. for /3 = 0.70 and L = 62. The coefficient b is related to the 
tunneling probability of the system between different topological sectors. Since most of 
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Q 


Q 




FIG. 2. Top: The free energy of the smeared topologieal eharge F {Q^) estimated by metadynamics 
and its statistical uncertainty, shown as a (orange) band at (3 = 0.70 with L = 62. The average is 
performed either over 30k sweeps (top left) or over 300k sweeps (top right). Bottom: comparison 
between the free energy F [Q^) estimated by metadynamics and —log{P{Q)) estimated by the 
standard HMC. The error on HMC results is not shown. 
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FIG. 3. (a) and (b) Time autocorrelation functions in logarithmic x-scale for the magnetization 
and the topological charge. The results correspond to jd = 0.75 and L = 60, for which E = 
0.6872601(25), fc = 5.160(3) and L/ic ~ 11-63. 




FIG. 4. (a) and (b) Topological charge as a function of the number of sweeps with the standard 
HMC and metadynamics respectively at jd = 0.75 and L = 60. In the calculation of v we have 
conventionally eliminated jumps in which the system returns to the original value of the topological 
charge in less than 10 sweeps. 


the transitions occnr with AQ = ±1, we expect that b ~ Kexp[S'/], where k is an entropy 
factor which increases with the volnnie and Sj is the one-instanton action Sj = 2 ttN fd. 
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FIG. 5. Scaling of the frequency v as a function of with metadynamics (dashed) and HMC(full 
line) for three different values of the case of the standard HMC, for af,G ^ 9 the changes 

of the topological sector are so rare that only an upper limit can he estimated. 

B. The central question: Towards the continuum limit 

In order to demonstrate that the approach presented in this work allows addressing 
efficiently the problem of topological trapping, we now discuss the scaling of the autocor¬ 
relation time as a function of the lattice spacing a ~ namely of ^g, and of the physical 
volume L/f^G- In Fig. and Tab. |IVB| we display the dependence of v on f^G with HMC and 
metadynamics at several values of As expected, because of entropy, in the standard 

HMC V is an increasing function of the lattice size, and this corresponds to an increase 
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in the dispersion of the topological charge, namely of (Q^). At larger values of (3, as we 
proceed toward the continuum limit, with the standard HMC, v decreases exponentially 
as a function of the correlation length EEI, so that for ^ 9 the changes of the 
topological sector are so rare that only an upper limit for u can be provided. 

The fact that the system is locked in a given topological sector may severely bias the lattice 
predictions for several important physical quantities, in primis the topological susceptibil¬ 
ity. The apparent small errors with HMC for the topological susceptibility at large values 
of (3, for example at (3 = 0.80, with L = 80, see Tab. are illusory because they refer 
to an average taken essentially at hxed topological charge and thus its value is affected 
by a systematic error. This is dramatically clear by comparing the values of the renor¬ 
malization group invariant product Xt^c column of Tab. A and [^ /? = 0.8 is 

just the largest value at which we could made a comparison between standard HMC and 
metadynamics, since above that it is impossible to get sensible results for the topological 
susceptibility using HMC. Indeed, with the standard HMC the number of sweeps necessary 
to get the correct results increases exponentially with (3 while in metadynamics increases 
only linearly. 


The behavior observed with the HMC is to be contrasted with the results obtained with 


metadynamics (corrected for the bias using eq. 25) since in this case u is sensibly flatter, 
and the simulation spans all the possible sectors of Q allowing to produce reliable results. 


We hrst analyze the scaling of the performances with the volume. It is known that in 
HMC the autocorrelation time of topological charge does not increase with the volume, 
as the increased tunneling rate (due to larger entropy) compensates the larger range of 
topological charges to explore. In metadynamics it is a priori not clear how the mechanism 
of tunnelling is affected by the volume. We note empirically (see Fig. that the tunnelling 
rate does not change with the volume, so one would expect that the algorithm efficiency 
degrades (mildly) with the volume. This is partially sustained from the observation of 
the errors reported on the last two columns of Tab. Nonetheless we remark that, no 
matter which volume is chosen, we still expect metadynamics to scale better with the 
lattice spacing than HMC. This means that for any physical volume there is a critical 
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lattice spacing below which using metadynamics is convenient. 


The hnal issue that we have to address is the scaling of the error of observable quantities 
measured in the metadynamics simulations. As we proceed to the continuum limit, the 
weight of conhgurations with non integer decreases, but the system keeps exploring 
all values G with equal probability after that the bias potential has 

converged. The unbiased expectation value is recovered by the reweighting procedure of 


eq. 25 The impact of the reweighting can be estimated modelling the scaling of the bias 


potential with the lattice spacing. Assuming that in the parametrization of eq. 26 the 
coefficient b, corresponding to the height of the barriers, grows approximately as 1/a, one 
expects a mild increase of the error due to reweighting. By comparing Tabs. [A| and [Bjone 
is comforted by observing that at small f3 (where HMC simulations are reliable) the errors 
for E and Xm obtained with metadynamics are comparable with those obtained with 
HMC, and increase very slowly with /3, not faster than the standard HMC. 


For the case of the renormalization invariant combination Xt^h observe that in the 
case of HMC the errors growth relatively less than with metadynamics, but the value drops 
dramatically at large due to the freezing of the topological charge. With metadynamics 
we hnd similar values across the whole range of explored lattice spacings, signaling that 
the algorithm samples correctly the distribution of the topological charge, as opposed to 
HMC. 


In spite of the great effort undertaken in the past |21 El E] even for the HMC algorithm 
it is not clear yet whether the topological modes are affected by exponential or power¬ 
like slowing down. In this hrst work we do not aim at determining accurately the scaling 
properties of our current implementation of metadynamics. We can however note that 
the observed growth of the error on the topological susceptibility as a function of the 
lattice spacing is compatible with a power law divergence of topological autocorrelation 
time, with an exponent that we estimate to be 2 -G 3. This has to be compared with the 
pure HMC simulations, for which the exponent of the power-like ansatz has steadily been 
suggested to be close to 5. We consider this a substantial improvement, that allows already 
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to have reliable results at (3 much larger that the HMC. Based on previous experience |33) . 
further improvement might be obtained including a larger set of collective variables to the 
algorithm, possibly coupling including in the bias other slow modes. 
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~ qN 

Ns 

L 

Cg 

Life 

E 

Xm X 10^ 

Xt X 10^ 

Xtfl X 10^ 


3M 

32 

2.7009(13) 

11.848(15) 

0.799324(9) 

12.182(3) 

11.66(18) 

85.1(1.3) 

0.65 

3M 

46 

2.7024(19) 

17.022(12) 

0.799352(8) 

12.148(3) 

11.91(26) 

87(2) 


3M 

64 

2.7075(23) 

23.64(2) 

0.799357(4) 

12.147(2) 

11.6(3) 

84.9(2.2) 


3M 

44 

3.748(3) 

11.739(8) 

0.739053(4) 

19.613(8) 

5.56(24) 

78(3) 

0.70 

3M 

62 

3.742(3) 

16.57(5) 

0.739072(4) 

19.545(4) 

5.93(14) 

83(2) 


3M 

88 

3.743(4) 

23.508(27) 

0.739066(3) 

19.543(5) 

5.88(26) 

82(4) 


3M 

60 

5.169(7) 

11.608(16) 

0.687256(6) 

31.890(31) 

2.5(3) 

67(8) 

0.75 

3M 

84 

5.154(8) 

16.299(24) 

0.687270(6) 

31.725(21) 

3.7(5) 

97(13) 


3M 

120 

5.152(16) 

23.28(7) 

0.6872678(21) 

31.723(18) 

2.9(5) 

78(13) 


3M 

80 

7.141(13) 

11.203(20) 

0.642261(4) 

52.48(6) 

0.74(12) 

38(6) 

0.80 

3M 

114 

7.102(8) 

16.051(18) 

0.6422752(22) 

52.120(23) 

0.29(6) 

15(3) 


3M 

160 

7.075(20) 

22.61(6) 

0.6422784(15) 

52.044(22) 

0.25(8) 

13(4) 


TABLE A: Parameters and measured quantities of the stan¬ 


dard HMC simulation. All the quantities are measured every 
5 sweeps exeept the energy and the topologieal eharge whieh 
are measured every sweep. 


~ qN 

Ns 

L 

^G 

L/^g 

E 

Xm X 10"^ 

Xt X 10^ 

Xifh X 10^ 


3M 

32 

2.7000(20) 

11.852(9) 

0.799310(14) 

12.181(4) 

11.44(14) 

83.4(1.1) 

0.65 

3M 

46 

2.6994(24) 

17.041(15) 

0.799366(6) 

12.1448(29) 

11.83(19) 

86.2(1.4) 


3M 

64 

2.702(4) 

23.696(36) 

0.799354(6) 

12.146(3) 

11.32(15) 

82.6(1.1) 


3M 

44 

3.743(4) 

11.756(12) 

0.739049(10) 

19.608(10) 

5.68(11) 

79.6(1.5) 

0.70 

3M 

62 

3.742(6) 

16.576(16) 

0.739069(6) 

19.550(7) 

5.66(12) 

79.2(1.7) 


3M 

88 

3.750(5) 

23.47(3) 

0.739071(4) 

19.539(4) 

5.68(25) 

80(3) 


3M 

60 

5.149(7) 

11.653(15) 

0.687281(7) 

31.795(21) 

3.01(7) 

79.8(1.8) 

0.75 

3M 

84 

5.141(7) 

16.339(23) 

0.687270(5) 

31.717(13) 

2.91(14) 

77(3) 


3M 

120 

5.183(11) 

23.15(5) 

0.687272(4) 

31.735(10) 

3.2(3) 

86(8) 


3M 

80 

7.079(7) 

11.302(11) 

0.642289(5) 

52.141(30) 

1.78(19) 

87(10) 

0.80 

3M 

114 

7.098(19) 

16.06(4) 

0.642288(4) 

51.98(4) 

1.80(23) 

89(12) 


3M 

160 

7.072(19) 

22.62(6) 

0.6422797(26) 

51.971(22) 

1.79(20) 

88(10) 

0.85 

3M 

108 

9.672(18) 

11.166(21) 

0.602879(5) 

85.91(8) 

1.00(16) 

91(15) 

0.875 

3M 

126 

11.073(21) 

11.379(21) 

0.584949(4) 

110.72(14) 

0.69(16) 

77(20) 

0.9 

3M 

148 

13.32(4) 

11.112(30) 

0.5680742(29) 

143.28(22) 

0.42(23) 

73(41) 


TABLE B: Parameters and measured quantities of the meta- 


dynamies simulation. 
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/? 

L 

1/ HMC 

V metadynamics 

0.65 

12 

(6.4 ±0.1) X 10"3 

(11.5 ±0.4) X 10-3 

0.70 

12 

(7.7 ±0.8) X 10-^ 

(5.9 ±0.2) X 10-3 

0.75 

12 

(5.0 ±0.3) X 10-5 

(2.2 ±0.1) X 10-3 

0.80 

12 

(1.8 ±0.6) X 10-5 

(8.7 ±0.5) X 10-^ 

0.83 

12 

(2.4 ± 1.4) X 10-^ 

(6.2 ±0.2) X 10-^ 

0.65 

18 

(11.1 ±0.1) X 10-3 

(13.4 ±0.2) X 10-3 

0.70 

18 

(1.1 ±0.1) X 10-3 

(5.2 ±0.1) X 10-3 

0.75 

18 

(9.3 ±0.5) X 10-5 

(1.84 ±0.04) X 10-3 

0.80 

18 

(2.5 ±0.9) X 10-5 

(8.2 ±0.3) X 10-“^ 

0.65 

25 

(16.2 ±0.2) X 10-3 

(16.4 ±0.2) X 10-3 

0.70 

25 

(2.9 ±0.2) X 10-3 

(4.62 ±0.06) X 10-3 

0.75 

25 

(1.7 ±0.1) X 10-^ 

(1.30 ±0.04) X 10-3 

0.80 

25 

(7.5 ±4.3) X 10-5 

(7.1 ±0.3) X 10-^ 


TABLE C. Values of the frequency u as a function of the coupling for different physical volumes. 
In the table we only give the values of v that have been reliably determined, excluding those corre¬ 
sponding to f3 so large that, in the case of HMC, only an upper bound can be given. 


V. CONCLUSIONS 


In this paper we have shown that the metadynamics approach [261 EH can be used 
to simulate CP^~^ improving dramatically the problem of the slowing down observed in 
numerical simulations for quantities related to the topological charge. In particular we 
have studied the iV = 21 case, showing that we are able to reconstruct the free energy of 
the topological charge, F{Q), as a function of the coupling and of the volume. 

The much reduced slowing down allows us to study a range of fi much larger than that 
available with ordinary HMC. Further improvement might be obtained by biasing a larger 
set of collective variables. 
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It seems straightforward to extend the general method exposed in this paper to the case 
of QCD. It will be interesting to investigate whether metadynamics will also work in this 
case. 


The possibility of measuring of F{Q) suggests that we can further improve the efficiency 
of the algorithm by using the information about its form as obtained in a preliminary run, 


cfr. eq. (26). This would allow simulating our system by using as importance sampling the 
knowledge of free energy itself added to the original action of our theory. 


The knowledge of the free energy of the topological charge {Q^ in our case) is very 
important for the study of relevant physical quantities. Actually, if we know the free energy 
of the topological charge, then we can compute the expectation value of any physical 
quantity at arbitrary values of ^-vacuum, and not only close to the origin, by simply 
averaging them with the appropriate weight ~ cos(6'Q) where we used the property 

that F{—Q) = F{Q). Metadynamics might thus also offer a solution to the difficulties 
encountered in simulating theories with complex actions (for example QCD with a non 
zero 6 term or at hnite chemical potential). 


Appendix A 


Taking inspiration from QCD we dehne recursively the (n + l)-level stout smeared link 
in terms of the n-level stout smeared link A” as 


'm;// 


\(n+l) _ \(n) 

'^m:u ^ '^m\u ’ 


'm;fi 


m]fi 


where 


and Sm]^i are the staples 


n. 


= Tm I 




c(n) _ [xH \(n) ,(n)* .(n)* .(n) .(n) 






26 


being p a small real number. The recursive procedure is based on the 0-level stout smeared 
links, that are nothing but the original (non-smeared) ones . 

As an effect of the exponentiation, the level-1 link results to be the average of the 
original link and the collection of all the infinite paths surrounding the two nearby 1x1 
plaquettes: 


+ 


P 




The definition is recursive and the n-level smeared link will involve contribution of links 
as far as d < n sites (each of them damped by a power Averaging a link with paths 
built in terms of its neighbors suppresses ultraviolet fluctuations, and reduces the noise in 
observables built in terms of the links. The parameters p and n can be tuned separately, 
but what indeed matters is the product pn. We found it convenient to use n = 2 and 
vary linearly p as a function of (3 {p = 0.08 at (3 = 0.65 and p = 0.13 at (3 = 0.90) which 
allows to separate the distribution probabilities of for different topological sectors as 


explained in Sec. Ill 
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